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Numerical Techniques in Li'.ear Duct 
Acoustics - 1980-81 Update 
Kenneth J. Batimelster 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 

A review is presented covering finite element and finite difference 
analysis of small amplitude (linear) sound propagation In straight and 
variable area ducts. This review stresses the new work performed during 
the 1980-1981 time frame, although a brief discussion of earlier work Is 
also Included. Emphasis Is placed on the latest state of the art In 
numerical techniques. 

NOMENCLATURE 

Cp ambient speed of sound, m/s 

■* 

c^ group velocity, m/s 

★ 

dp duct diameter, m 

F Initial condition vector, equation (22) 

f* frequency, Hz 

t function of y, see equation (18) 

H* duct height, m 

I number of axial qrld points 

. fTT 

•• number of transverse grid points 

K coefficient matrix, equation (22) 

L* length of duct, m 

Nju(j number of storage locations in sub matrix 


n 


normal 
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pressure, N/m 

* *2 

time dependent acoustic pressure, 

spatially dependent acoustic pressure, equation (3) 

wave envelope pressure, equation (25) 

radius of duct, m 

dimensionless time, t*/tp 

period 1/f*, sec 

time step 

time dependent axial dimensionless acoustic velocity, 

U(x,y,t). U*/c* 

Mach number, Uq/Cq 

spatial dimensionless axial acoustic velocity u(x,y), u*/Cq 
time dependent dimensionless transverse or radial dimensionless 
acoustic velocity, V(x,y,t), V*/Cq 
spatial dimensional transverse or radial dimensionles acoustic 

velocity, v*/c^j 

axial coordinate, x*/H* or x*/r* 
axial cirid spacing 

dimensionless transverse coordinate, y*/H* 
transverse grid spacing 
impedance, kg/m‘ sec 

dimensionless specific acoustic impedance 
dimensionless frequency, H*/f* c^ (cartesian) or r^/f* c^ 
(cylindrical) 

dimensionless specific resistance 
dimensionless wave length 
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ambieni air density, kg/m^ ^ 

dimensionless acoustic potential function, 
dlmenslonless reactance 
angular freauency (?»f*) 

Subscripts: 

e exit condition 

1 axial Index (fig. 1) 

j transverse Index (fig. 1) 

0 ambient condition or mean flow variables 

* dimensional quantity 

k time step Index 

(1) real part 

( 2 ) imaginary part 

INTRODUCTION 

"Steady'' state finite element theories and transient finite difference 
theories are currently being applied to the design of loud speakers, ven- 
tlllatlcn systems, mufflers, and the acoustically treated nacelles of turbo- 
jet engines. In this paper, a review Is presented covering both the finite 
difference and finite element analysis of small amplltuoe (linear) sound 
propagation In straight and variable area ducts with and without flow. 
Reference 1 contains a very recent literature review of the techniques, 
advantages, limitations and applications associated with the various numeri- 
cal solutions of the sound propagation equations In ducts. In an extension 
of reference 1, this review emphasizes the new work performed during the 
1980-1981 time period and the state of the art of numerical simulation of 
sound propagation In ducts. 
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Application of numerical theories to room acoustics, structural acoustic 
radiation systems and tne determination of duct eigenvalues will not be 
included in the detailed discussions. However, since these research areas 
are often of interest to acousticians studying duct acoustics, these works 
have been included in the list of references (refs. 2-22) for completeness. 

The numerical theories to be considered in detail in this paper employ 
full two or three dimensional finite element or difference tneories which 
have provisions for inlet and outlet flow such as in a muffler or jet engine 
duct. To delineate the starting point of the present upaate, the literature 
cited in reference 1 is again tabulated herein (refs. 23 to 63). These 
references will also be cited in conjunction with brief summaries of earlier 
work. References 64 to 86 represent the new literature introduced during 
the 1980 to 1981 update period. 

First some background information is briefly discussed relating to the 
properties of sound propagation in ducts and how they impact on numerical 
solutions. Next, the various numerical works developed in the 1980-1981 
period are presented in the following order: 

(1) Steady State Finite Difference Theory 

(2) Steady State Finite Element Theory 

vi) Special Numerical Transformations 

(4) Transient Numerical Theory 

BACKGROUND INFORMATION 

In preparation for the new (1980-1981) literature on finite oifference 
and finite element techniques as applied to duct acoustics, the simplest 
versions of the sound propagation equations and boundary conoitions a-'e 
first briefly revieweo and some inherent problems associated with numerical 
solutions of the sound propagation equations are high-lighteo. 
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(1/n may also be thought of as a dimensionless speed of sound, in this 
case, equation (1) corresponds more closely to the dimensional form of 
the wave equation). The starred quantities represent dimensional vari- 
ables. All symbols are defined in the nomenclature. The sound source is 
generally assumed to be harmonic in nature and to vary as e (in 
dimensionless form as In this case, the pressure and acoustic 

velocities will all be functions of If the acoustic pressure 

is assumed to be of the form 

P(z.y,t) = p(x,y)e’^’''^ (3) 

the wave equation (1) becomes 

?!p + ifp + (2,n)^ p = 0 (4) 

which is the classical Helmholtz equation. 

Solutions to equation (4) are called "steady" state solutions because 
the equation is independent of time. In this case, p(x,y) actually re- 
presents the Fourier transform of P. On the other hand, solutions to 
equation (1) are time dependent or transient. Once the initial start up 
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condition has died, tne transient and “steady" state solutions are simply 

related by equation (3). 

Governing Equations with Flow 

The wave equation (1) takes on slightly different forms tor tne cases 
of uniform mean flow ana for irrotational flow (see ref. 40 ). For uniform 
mean flow the wave equation becomes 
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(b) 


As is equation (1), this form of the wave equation is ideally suited for 
numerical analysis in that only one dependent variable is involved in tne 
solution of the propagation equations. The wave equations are obtained 
by combining the linearized continuity and momentum equations which 
describe the propagation of small acoustic disturbances. Unfortunately, 
for realistic shear flows, the continuity and momentum equations cannot 
be combined into a single second order wave equation. In these cases, 
the continuity and momentum equations must be solved simultaneously. For 
parallel shear flow in a two-dimensional rectangular coordinate system, 
these equations are 

3P.__LliL.l2l-!i2.iL (6) 

3 t * ~ ^ 3n ” n ay n ax 

at * ~ n ax ~ n ax n ay 

- - i if. - — — 

at “ ” n ay ” n ax 


For “steady" state, these equations become 


( i2lTTl)p = 
(i2irn)u = 


9U 

3V 

- u 

i£ 

' * ' (9) 

3X * 

ay 

0 

3X 
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(10) 
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ay 
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(11) 

ay 

0 

ax 




As will be shown later in detail, the "steady" state solution of the 
sounc propagation equations requires storage of large matrices to obtain 
a solution. Consequently, the inclusion of three dependent variables in 
the shear flow problem requires an order of magnitude increase in the 
storage capability and solution times of a computer over a similar prob- 
lem in the absence of shear. Therefore, to obtain optimum storage and 
run times, separate programs are required for no flow, irrotational 
flows, and shear flows. Fortunately, as will be shown shortly, the 
transient solution to the governing equations does not require the stor- 
age cf large matrices. Therefore, the generation of a single optimum 
program to cover the complete range of no flow and shear flow seems a 

reasonable goal. 

Discretize the Continuum 

In either the finite difference or finite element numerical analysis, 
the continuous acoustic flow field is lumped into a series of grid points 
or elements as shown in figure 1. Instead of obtaining a continuous 
solution for the acoustic pressure, for example, its values are c^btained 
only at isolated node points. In the finite difference analysis, a sinv- 
ple rectangular pattern of discrete points is almost always used. The 
finite difference approach is generally restricted to uniform ducts 


- - 
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unless special bookkeeping or transformations are employed. On the other 
hand, a great variety of finite element patterns (see table I) can be 
easily employed to handle complex geometric variations in duct walls. 

Two such element patterns are displayed in figure 1. The more nodes per 
element the fewer will be the number of elements needed to resolve the 
acoustic Tield. However, since the sire of the solution (global) matrix 
is proportional to the number of nodes, considerable more computational 
time and computer core memory are required for the higher order elements. 
Grid Point Requirements 

in sound propasatlon In ducts, tha acoustic pressures and velocities 
oscillate do«n the duct. The puestion naturally arises as to how oany 
prid points or elements are reouired to accurately resolve these oscilla- 
tions. for example, consider a hard wall duct, semi-infinite in extent 
with a sinusodtal plane wave pressure at x . o, and a uniform Hath num- 
ber U„ in the duct. The one dimensional -steady state plane wave 
solution of euuation (5) in conjunction with efluation (3) yields 

i2irnx 

(1) ^ t J2) _ 

P s= P ^ » P - ^ 

Thus, the real and imaginary components of the acoustic pressure oscil- 
late in the axial direction with a wave length 


The number of axial grid points to obtain accurate acoustic pressure 
and velocity profiles may be estimated from the following -rule of thumb' 

given in references 24 and 25 



I 


X 


(14) 


Thus, for the unit frequency (n « 1), unit length (L*/H* . 1) ana no mean , 
flow, 12 grid points are necessary to describe adequately tne sinusoidal 
form of the spatial pressure dependence, as shown in figure 2. If the 
frequency or length is doubleo, or the Mach number is -0.5 (flow 
towards the source], clearly twice as many grid points will be required 
to describe the wave, since two wave lengths of sound must now be re- 
solved. Vihen higher order transverse acoustic modes are present, the 
axial wave length of sound will be longer and the required number of 
axial grid points will oe significantly reoucea. Therefore, equation 
(13) represents a conservative estimate of the grid point requirements. 

Tne number of transverse grid points will oepehd on both the source 
transverse modal distribution and the generation of higher order modes in 
the duct. To resolve ail the higher order propagating mooes that could 
occur in tne most general case, the number of grio points in tne trans- 
verse direction would oe proportional to the frequency n (rets. 61 and 
63). The total number of grid points would be the product of the points 
\ in the axial and transverse directions. In high frequency (large tj) 

cases, the required number of grid points or elements can become very 
large. 

In determining the required number of finite elements, however, tne 
constant 12 in equation (6) will oe smaller because the interpolation 
function in finite element theory is generally of higher order than the 
first order linear difference approximation used to establish equation 
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(7). For a particular finite elenent, some numerical -experiilientation 
will be rcquii'ed to determine it sufficient numbers of elements have been* 
employed. 

Wall Boundary Conditions 

The boundary condition at a hard wall duct is simply that the normal 
acoustic pressure gradient is zero 


s IE 
an “ an 


0 


(15) 


or that the normal acoustic velocity is zero. 

- 0 (16) 

For a soft wall (absorbing) duct, an impedance is specified at the duct 
wall 

c » 4^- ix (17) 

fo’^o 



> . 






St.- 


The wall resistance is represented by e^. while the reactance by 
X. Both and x a**® assumed known input quantities. 

In the "steady" state analyses, application of equation (17) is 
sTraioht-fcrward and offers no computational problem. In the transient 
analysis, however, direct application of equation (17) can lead to 
numerical instabilities. The integration technique presented in refer- 
ence 61 seems to alleviate this problem. 

Entrance Conditions 

For a rectangular duct, the boundary condition at a liner entrance 
P(o.y.t) can be of the form 

P{o,y.t) - f{y)’“*^* (18) 
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where f(y) represents the transverse variation in acoustic pressure. 

Similar eapressiona tho vclacity pptentia! or acoostto volocit. cap be . 
used. Id cylindrical coordinates, the radial coordinate r «oold replace 

, i„e,uation(l8). Often the function form f(y) is associated with the 
normal-mode analytical solution-cosines for the rectangular coordinate sys- 
tem and Bessel functions tor the cylindrical coordinate system. 

Physically, boundary. condition (18) represents the sum of a foruard and 
reflected acoustic »ave at x . 0. In a uninuely different approach, 
Eversman, et al. (ref. a?l, assumed a uniform infinite hard wall duct 
upstream of the duct section of interest and that eouation (18) represents 
only the forward propagating wave. The amplitudes of a truncated scries of 
normal mode reflection coefficients were determined by matching the pres- 
sure and pressure gradients in the infinite duct to the finite element 
ualues at the entrance. This form of termination is particularly useful m 
conjunction with many experiments which employ an anechoic entrance or exit 

ccndition. 

Exit Condition 

In applying the finite difference or finite element analyses to a 
turbojet enolne inlet, for example, the grid system has generally been con- 
fined to the internal portion of the inlet. Thus, the engine has been 
modeled as a short duct. In this case, an impedance of the form 

p(L*' (15) 

S = VTTTT 


is often used to close the boundary, 
approximate an anechoic termination, 
ance for single modes in an infinite 


The value of Cg 1= chosen to 
General values for the exit imped- 
hard wall duct (refs. 38 and 40) are 


11 


often used for the anechoic approximation, for arbitrary multi-modal' 
wave forms, however, a geheraV impedance eqeaticn is not available for an 
exact simulation of an anechoic termination. 

As with the entrance condition, Eversman, et a1. (ref. 22} use the 
tecnnique of modal decomposition to obtain an accurate simulation of an 
anechoic termination. 

Reference 24 (Appendix E) suggests another possibility for simulating 
a non-reflecting interface at the duct exit. By aoding on an additional 
length of duct with acoustical damping, the reflections from the duct 
exit iSfill oe effectively damped before they can re-enter the orininal 
portion of the duct. This technique has worked extremely well in soft 
wall duct problems (ref. 24) and is presently being employed in the 
transient analysis. Also, in the transient analysis, an anechoic termi- 
nation can be modeled exactly by extending the duct such that the steady 
state solution is Obtained before the arrival of a reflected wave from 
the duct exit. This approach, nowever, is costly from the standpoint of 
CGP:'ijter storage and run times. 

All the previous cases attempted to eliminate or at least reduce 
ref Uctions at the duct exit. However, in many practical applications, 
such as a turbofan inlet, reflections could be important for certain 
mooes. Consequently, continuing the gr'd structure from inside the duct 
into the far field (ref. 31) would simulate the actual dynamic process 
occurring at the lip of the duct. In the far field, all ^.jct mooes 

Hr * 

propagate and nave identical OqC^ exit impedance far from the 

exit. Tneretore, the closure problem is simplified out at the expense or 

a greater number of elements and a larger solution matrix. To eliminate 
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this l?rge matrix requirement, Kwgawa, et ai. (ref. 4iO) developed a conv- 
Dination of finite element and analytical methods to analyze -§ound prop.a- 
gaticn into the far field. 

Methods of Solution 

The methods of solution for the various "steady" ano transient 
theories will be discussed in detail later. Generally, the "steady" 
state finite difference and finite element numerical algorithms have been 
limited to low frequencies (small n) because of the grid point require- 
ments and more importantly because of tne large matrices associates with 
the solution of tne time independent equations. Special transformations 
have also been developed (refs. 54-60) to either reduce tne size of the 
"steady" state matrix or to eliminate it. As an alternative to the 
"steady" state theories, time dependent numerical algorithms (refs. 
61-63) were developed which eliminate tne need to store a large solution 


matrix. 


STEADY STATE FINITE DIFFERENCE THEORY 


In tne finite difference theory, the governing equation is written in 
difference form at each grid point. For the previous no flow example, 
equation (3) can be written in difference form by the usual 5-point dif- 
ference approximations. 



(20) 


or for AX equal to ay 


Pi-l.J * "i J-1 "i,j ^ '’i.rl ° 


( 21 ) 
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Equation (21) applies away from the Cell boundaries. For grid points on 
the boundary, the hard {eq. (15)) or soft wall (absorbing) impedance 
condition (eq. (17)) is imposed. Methods for generating the difference 
equations on the boundary are fully documented in references 23-28. 

Matrix Solution 

Tho collection of the various difference equations at each grid point 
forms a set of simultaneous equations that can be expressed as 

[Kl (PI - {F} (22) 

where [K] is the known coefficient matrix, (|T} is the pressure vector 
containing the P. . unknowns, and (T) is the known column vector 
containing the various boundary conditions. This matrix is complex 
becaust? of the complex nature of the source and impedance boundary condi- 
tions (eq. 17). As seen in equation (21), the frequency term (wnAx)^ 
subtracts from the main diagonal element such that the coefficient matrix 
is not diagonally dominant. As a result, iteration solutions cannot gen- 
erally be used. Equation (22) is usually solved by some elimination 
technique. Unfortunately, the elimination technique requires the storage 
»'f all the matrix elements. 

The total storage will be proportional to the square of the total 
number of grid points. Because of the banded nature of the difference 
matrix, sparse matrix techniques have been employed to reduce computer 
storaae and run times as much as possible. Quinn (ref. 28) partitioned 
the matrix K into tridiagonal form, which reduced the elements in the 
submatrix storage to 
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(23) 


'^sub " ■■ 

which led to a total storaoe requirement of 


'mpiex 


3456 n'^(L*/H*) 


(2<l) 


For n ■ 10, L*/H* « 1 and ■ 0, the total storage will be 
3.456x10^. which is quite_lariie. Higher frequencies and flows will 
greatly enlarge the storage. These submatrices are read in one at a time 
to obtain a solution of the matrix. Thus, the input-output time of the 
computer solution can also become very large. 

Update 1980-1981 

Finite difference solutions (ref. 64) were used to study plane wave 
sound propagation in a curved bend joined to two straight sections ot 
duct. A rectangular mesh and a cylindrical mesh were used in the 
straight and curved sections of the duct respectively. At the dis- 
continuities in the coordinate systems, a linear interpolation was used 
in the radial direction in conjunction with unequal meshes in the axial 
direction to derive the appropriate finite difference approximation. 
Eouation (Id) was employed at the entrance condition while moda: de- 
composition and matching were used for the anechoic exit condition. 
Unfortunately, the algorithms for the difference methods are not pre- 
sented in reference 64. For a variety of curved ducts, the numerical 
results aoreed closely with the experimental data for the reflection 
coefficient and the pressure amplitude across the duct. For frequencies 
just below the first cut-on mode frequency, the magnitude ot the cut-off 
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modes generated at tne interface were found to be greater in magnitude 
than the incident plane wave. 

With the c;:cepticn of reference £4, fcasic research on the “steady* 
state finite difference method appears on the wane. The bulk of recent 
"steady" state acoustic analyses employ the finite element theory. The 
finite element theory can readily handle duct area variation or curvature 
changes without resorting to special bookkeeping procedures. However, 
finite difference theory isagaln receiving renewed interest In the 
transient solution of the wave equation. This topic will be covered 
later in the paper. 

Status 

A finite difference program Is available In the literature for han- 
dling rectangular and circular ducts with uniform flow. Reference 28 
contains a complete list of fortran statements plus example problems to 
illustrate tne use of the program. Tne program also has provisions for 
analyzing variable area ducts if a mapping function is available. The 
mapping function for a cone and hyperbolic horn are listed. As written, 
the program is limited to 100 axial grid points and 20 grid points in tne 
transverse directions. This limited number of transverse gria points 
restricts the number of higher order modes which can be resolved. There- 
fore, some modification or the program would be necessary in the analysis 
of high frequency sound. 

In general for ducts with variable area, the finite element theory is 
tne most convenient to use. However, as pointed out by Quinn (ref. 53), 
the aavantage of finite elements is not as great when a mapping function 
is used to compute the mean flow field. In such a case, where the con- 


16 


formal transformation Is availaOie, the use of finite ouference tneoi*y 
would have an advantage over finite element theory. In this case, the 
program cited in reference 28 could be used very effectively. 

"STEADY" STATE FINITE ELEMENT THEORY 

Finite element theory is reviewed in depth in many references 29-50 
as well as the previous review paper (ref. 1). Consequently, only a very 
brief introduction to the finite element theory is now presentea. 

In the finite element theory, an appropriate element with its asso- 
ciated nodes is distributed in the duct, such as shown in figure 1. At 
each node labeled i, an unknown value of pressure p^ is assigned. 

The acoustic velocities would also be assigned to the nodes in the more 
general shear flow problem. Table I shows some additional elements cur- 
rently being used in acoustic studies. In contrast to the finite differ- 
ence theory, interpolation functions are used to determine the value of 
pressure between the nodes at any position inside the element. Again, 
taoli 1 lists some typical interpolation functions. 

Generally, Lagrange polynomials are used for interpolation when only 
the magnitudes of the dependent variables, such as pressure, are desired 
at eacn node. This is commonly called the C* continuity problem. In 
inlet ducts, Astley ana Eversman (ref. 51) showed that C* continuity 
leaos- to smooth solutions for the acoustic pressure proviaed a suffi- 
ciently small mesh is employed. If the mesh size is fixed, however, ana 
the Macn numoer is increased, the solutions for pressure become cusped 
(discontinuous) at each node. Improvements in resolution for the high 
Mach number cases require either greater mesh refinement with a corre- 
sponding increase in dimensionality or the introauction of more sophisti- 
cated Hermitian elements. 
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When using Hermite polynomials, items 1 and 6 in table I, the values 
of slope are also calculated at each node, commonly called the C prob- 
lem. In this case, f^wer elements are required to give smooth solutions 
at the nodes. However, C* continuity is more difficult to construct, so 
with the exception of items 1 and 6 in Table 1, Lagrange polynominals are 
used. Also, C’ increases the unknowns by a factor of 4 (item 6, table 
I), and the resulting solution matrix will increase by about 16. 

Next, either the Galerkin, least squares or variational finite ele- 
ment methods is used to constrain the pressures at the node points such 
that they satisfy the Helirioltz equation (3) or a more general form of 
the wave propagation equation. Generally, tne variational techniques 
have been used for problems without mean flow, while the Galerkin formu- 
lation is used when mean flow is present, see table I. Then the ele- 
mental equations at each node are combined into a general matrix (global) 
similar to equation (22). Again, this matrix must be solved to determine 
the unknown pressure p^. 

Matrix Solution 

With the exception of the least squares approach, the global matrix 
must be solved by some elimination technique just like the finite differ- 
ence solution. For the larger problems, an out-of-core banded solver is 
generally required with a moderate amount of in-core storage but much 
more input/output time. To reduce computer storage and run times, the 
nodes should have been labeled in a manner to obtain the minimum possible 
banded matrix. 
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Update 1980-1981 


No Flow 

In reference 65« Shepherd and Cabelll presented numerical and experi- 
mental techniques to determine the acoustic characteristics of dis- 
continuities in rigid curved wall duct systems when several modes propa- 
gate. • The finite element technique was used to solve the two-dimensional 
Helmholtz equation and determine the acoustic behavior of arbitrarily 
shaped duct elements in terms of modal amplitudes. Eignt-noae iso- 
parametric rectangular elements were used to describe the solution wUn 
haro wall boundary conditions. Shepherd and Cabell i applied tne method 
to a 90* mitred bend without mean flow and obtained good agreement 
between numerical and experimental results. In references 06 ana 67, 
Cabelli and Shepherd extended the experiments and theory of reference 63 
to a range of mitre bend geometries. The influence of change in geometry 
on the acoustic characteristics of 90* bends was established over a range 
or inner ana outer radii and the effects of including a turning vane were 
examined. 

Figure 3 displays a layout of the experimental apparatus used by 
Shepherd and Cabelli (ref. 65) to measure tne propagation of higner oruer 
duct modes through mitred bands. In the finite element simulation or 
this experiment, tne modal match boundary •conditions developed by Astiey 
and Evorsman (rei, 50), inrinite duct modal solutions coupled to finite 
element solution, were used to model the anechoic entrance ana exit 
condition. Good agreement between experiment aiio tneory was Obtained. 

In contrast, it would be difficult to employ an exit impedance boundary 
condition, equation (19), in this problem because of the aifticulty or 
specifying an impedance when multiple modes are present. 


In addition to the theoretical studies, soiTie duct experiments were 
conducted which were compared to previously-developed finite element pro- 
grams. In reference 68 and its improved version reference 69, the ab- 
sorbing properties of a linear locally reacting honeycomb liner were 
measured for an incident plane wave (no mean flow). Tne 5.08 cm square 
anechoic test section supported plane wave propagation up to about 
3.5 KHz; therefore, the standard plane wave p*c* exit imped- 
ance was useo to terminate the solution. Tl:e impedance of tne liner was 
determined in a separate normal incidence measurement. Very good agree- 
ment between finite element theory and experiment was observed at 0.9 KHz 
while ^n the high frequency tests of 2.5 KHz, significant discrepancy 
existeo in noth the pressure level ano phase at the trailing edge of the 
1 iner. 

Because of the rotational nature of the blades in a turbofan jet 
engine, much of the emitted turbomachinery noise is concentrated in 
spinning acoustic modes. The NASA Langley Spinning Mode Synthesizer 
shown in figure 4 has the capability to generate individually dominant 
circumferential modal patterns. In reference 70, comparisons of tne 
transmission ano reflection coefficients with calculations based on 
finite element models (refs. 37 and 46) were made. In this case, the 
experiment was conducted without mean flow with hard walls and a dominant 
spinning mode of order (1,0). In general, reasonable agreement between 
experiment and tneory was obtained. 

Uniform Flow 

Ling and Hamilton (refs. 71 and 72) have developed a finite element 
model for steady uniform flow based on the "steady" state form of equa- 
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tion (b). Ihe finite element formulation is baseo on the Gaierkin methoo 
using cubic isoparametric elements. The model was also used to approxi- 
mate sound propagation in a converging-diverging duct with flow. Since 
the 3 Ug/ax term has been neglected in their analysis, a comparison 
with the more exact theories of references 45 to 50 will be required to 
establish the range of validity of their finite element program as 
applied to flow with area variations. 

Potential Flow 

Sigman and Zinn (ref. 73) have applied their potential flow finite 
element program (refs. 40 and 43) to the problem of describing wave 
propagation in axisymmetric nozzles. Agreement between values of nozzle 
admittances computed from finite element solutions experimental admit- 
tances are extremely good. 

Shear Flow 

Astley and Eversman (ref. 74) have extended and improved their 
earlier analysis (ref. 51) on the implementation of the finite element 
method for acoustic transmission through a non-uniform duct carrying a 
high speeo suDsonic compressiole sheared flow. A finite element scheme 
based on both the Gaierkin method and the resioual least squares method 
using an eight model isoparametric element was described. Multimodal 
propagation was investigated by coupling of the solution in the duct non- 
uniform section to modal expansions in uniform sections. As in refer- 
ences 48 and 53, they found that the residual least squares formulation 
performs poorly in comparison with the Galerken approach. At high Mach 
numbers, evidence was provided that multi-mooal interactions become 
important necessitating a proper handling of the radiation condition. 
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Astley, Walkington and Eversman (ref. 75) investigated tne computa- 
tional efficiency and possible improvements to existing finite element 
models. Numerical results were presented in an attempt to evaluate four 
computational schemes comprising all combinations of Galerkin and resid- 
ual least squares with both Lagrangian (C’) and Hermitian (C^) ele- 
ments. A computationally inexpensive one dimensinal model was used for 
most of the results. It was demonstrated that an inevitable consequence 
of the Galerkin schemes appears to be the presence of spurious oscilla- 
tions excited by insufficient number of elements or inaccurate matching 
at the boundaries. Also, the wave envelope technique, to be discussed in 
the next section, showed a significant increase in accuracy over tne con- 
ventional solutions when a fixed number of elements was employed. 

Abrahamson (ref. 76) attempted the development of a 3 dimensional 
finite element program to analyze acoustic propagation in an axisymmetric 
duct with circumferential variations in wall impedance. Galerkin or 
least squares element formulations comoined with Gaussian elimination, 
successive over-relaxation, or conjugate gradient solution algorithms 
were investigated. Unfortunately, all the techniques proved impractical 
for handling realistic three dimensional problems, because of both very 
large storage and run times. 

Status 

Nearly all the finite element investigations have been concerned with 
the theory and accuracy of various forms of tne finite element method for 
the calculation of the acoustic transmission in ducts. In most cases, 
only in-core implementation of the finite element metnoa has been em- 
ployeo. Abrahamson (ref. 46), however, has developed a large scale com- 
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putational scheme utilizing an out of core banded solver for handling ‘ 
sound propagation in ducts containing a compressible mean sheared flow. 

The program has the capacity to hanole relatively large frequencies 
(n “ 10), duct lengths, and Mach numbers. Both a computer program manual 
(ref. 77) and detailed flow charts (ref. 78) are available. The detailed 
Fortran listing can be obtained directly from the author. 

The program solves the separate linearized continuity and momentum 
equations using the Galerkin finite element method. As shown in fig- 
ure 5, a simple linear rectangular Lagrangian element is used for which 
the acoustic pressures and velocities are to be determined at each node. 
The numoering system snown in figure 5 leads to a global matrix wnich 
consists of 3 diagonal bands. For a typical two dimensional problem, the 
width of each bano is 9. The structure of the global matrix is the same 
for the uniform, non-uniform and variable area case shown in figure 5. 
These elements are stored in a packed matrix which ignores the large 
number of zeros in the general global matrix. The packed blocked tri- 
diagonal matrix is then solved in the usual manner by forward and back- 
ward substitution. 

At the present time, the program has been successfully applied to a 
variety of no-rlow and uniform flow problems for straight and variable 
area ducts. Additional checking and debugging may still oe necessary whe 
the program is applied in a sheared flow situation. It would be useful 
if the modal matching boundary condition could be incorporated into this 
program so that anechoic termination could oe conveniently simulated wnen 
higher order modes can propagate. 


'•STEADY*' STATE SPECIAL NUMERICAL TRANSFORMATIONS- 

In a turbofan inlet, high frequencies (n * 30 to 50) and flow Mach 
numbers lead to pressure and velocity oscillations where many axial wave 
lengths variations are present in the duct. The number of grid points or 
elements reauired to track these oscillations requires extensive computer 
storage and run times which prohibit the analysis of many duct problems. 
Also, in a computer optimization process, hundreds of calculations are 
often required to determine a desired liner configuration. Therefore, 
any reduction in the number of grid points or elements in a numerical 
analysis will significantly reduce the cost of obtaining the desired 
solution. As cited in reference 1, three special numerical techniques 
have been developed in order to reduce the storage requirement and in- 
crease the accuracy of the numerical solution. 

Wave Envelope Technique 

The wave envelope technique attempts to reduce the axial oscillatory 
part of the wave pressure profile by transforming the wave equation into 
a new form whose solution is non-oscillatory. For example, if one domi- 
nant mode is present in a duct, the pressure can be assumed to be of the 
form 

-i2itx 

P(x,y) = PQ(x.y) e ^ = P^Ix.y) (25) 

where p(x,y) is the oscillatory pressure shown by the solid line in fig- 
ure 6 (real part only) and p^ represents the pressure amplitude or 
envelope shown by the dashed line in figure 6. The amplitude in figure 6 
falls off because of absorption by the soft wall. Substituting equation 
(25) into the Helmholtz equation (4) yields the wave envelope equation 
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(26) 
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Solutions of equation (26) have lead to order of an magnitude reouc- 
tion In grid points and computer storage compared to conventional solu- 
tions of the Helmholtz equation, equation (4), as discussed In references 
54.47. The technique is limited to those cases where some reasonable 
estimate of the wavelength can be made a priori or by some slnpler analy- 
sis. Re-running a problem with Increased grid points will generally be 
necessary to check the validity of the solution. 

Marching Technique 

Although Initial value numerical solutions of elliptic equations are 
generally unstable, the two-almenslonal Helmholtz wave equation can oe 
solved as an Initial value problem using explicit marching techniques 
(refs. 58 and 59). Compared to standard finite-difference or finite ele- 
ment boundary value approaches, this numerical technique Is oraers of 
magnitude shorter In computational time and required computer storage. 
This technique Is limited, however, to high frequencies and to cases 
where reflections are small. 

Far Field Coupling 

To obtain far field acoustic radiation patterns and to eliminate the 
necessity of specifying the impedance of the duct at its exit, finite 
elements can be extended frcm the internal portion of a duct Into the far 
field (ref. 31, fig. 16). Generally, the increased computer storage 
requirements make this approach very costly. Kagawa, et al. (ref. 31) 
recently developed a combination of finir» ele«>ent (in ouct) ana analyti- 
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cal methods (Green-theorem-far field) to analyze sound propagation from 
an acoustic horn. Bauraelster and Majjigl (ref. 44) presented a matrix 
partitioning approximation that separates the uuct and tar field into two 
or more regions and thereby reduces total matrix storage. 

Update 1980-1981 
Wave Envelope Technique 

Astley and Eversnar. (ref. 79) have adapted the wave envelope approach 
to their finite element program which was documented in reference 74. In 
contrast to the earlier works which considered mainly plane wave applica- 
tions in uniform ducts, they considered multimodal propagation In soft 
walled straight and variable area ducts with flow. Astley and Eversman 
reported that significant Increases In accuracy result from a relatively 
crude wave envelope modification to conventional finite element schemes, 
in practical terms, as suggested In previous studies, Astley and Eversman 
agree that this technlpue could significantly reduce computer storage and 
tun times In realistic aeroacoustic configurations where the propagation 
of high frequency modes places great demands on the axial resolution of 
conventional numerical schemes. 

Far Field Coupling 

In oetermining far field radiation patterns from Inlets, generally 
the '‘Interior" region of the duct and unbounded "exterior" region of the 
inlet are treated separately. In the aosence of flow, Horowitz, Sigman 
and Zinn (ref. 80) developeo a hybrid iterative solution approach using 
finite elements for the interior of tne duct and the Integral technique 
for the exterior region in order to obtain a continuous numerical solu- 
tion for the acoustic field. In the first Iteration, an Impedance bouno- 
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ary condition was assumed at Jhe interface between the interior and 
exterior regions. The boundary condition was then used to obtain inte- 
rior and exterior (radfatibri) acoustic solutions. Since the prescribed 


boundary condition at the interface was an approximation, a discontinuity 

• f 

arises in the computed acoustic fields at the interface. In this case, 
both the interior and exterior acoustic fields were in error. This error 
was eliminated by an iterative matching of the computed acoustic fields 
at the interface which also provides the correct boundary condition at 
that location. 

Resul^ng solutions for several simple cases were compared with known 
exact- solutions and good agreement was demonstrated. In all cases, the 
convergence was very rapid. Under a NASA grant, work is continuing to 


adapt this matching procedure for the case when flow enters the duct. 
Status 

At the present time, these special techniques have not been in- 
corporated into a documented computer program. Therefore, if desired, 
individual programs must be revised to accommodate these special tech- 
niques just discussed. 

TRANSIENT FINITE DIFFERENCE THEORY 
The transient analysis begins with a harmonic (e " ) sound source, 
equation (18). at x » 0 radiating into an initially quiescent duct. 
Next, an explicit iteration solution of the wave equation, equation (1) 
is used to calculate both the transient as well as the "steady" state 
solution. The "steady" state is obtained when the value of 


P(x,y) = P(x,y,t)/e’^*" (27) 
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at the duct exit reaches a cohstant value. For comparison and reporting 
purposes, all the transient pressures are converted Into “steady" pres- 
sures according to equation (27) thereby eliminating the temporal 
variations. 

To Illustrate how the transient equations are programmed, the second 
derivatives in the wave equation, equation (1), can oe represented by the 
usual central differences in time (superscript k) and space (subscripts 
1, j) - - 

/„k+l .,pk ^ /p*^ 29^ + \ - 2P*^ • 



where at, ax, ay are the time and space mesh spacings, respectively. 

k+1 

Solving equation (28) for j yields 

pk+1 -pk pk-1 ^ / at + pl^ _ 4p*^ . + p*^ . , + p*^., .1 

J " \nax/L^1-l,j ^1,j-l i.J 


(29) 


where ax equals ay. The Iteration procedure Is explicit since all the past 

values of P^ are known as the new values of P are computed. Since 

equation (29) Is a simple algebraic equation, storage is only 

k k— 1 

required for tne solution vectors P^ j and ’s, matrix 

storage and manipulation are not required. 

As with any explicit iteration solution, the time step at must be cnosen 
sufficiently small to Insure numerical staollity. Nevertheless, tne transient 
solutions have been found to be an order of magnitude faster than tne similar 
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"steady" solution, ref. 61, for plane wave propagation. Using the von Neumann 
stability theory as a guide {refs.- 61-53), the -time dependent solutions have 
been found to be stable without flow and with uniform flow. When axial and 
transverse variations in the mean flow are considered, some numerical experi- 
mentation will no doubt be required to determine the maximum time increment 
associated with numerical stability. 


1981 Update ' ' 

As previously formulated in the literature (ref. 63), the transient method 
did not converge to the "steady state" solution for cut-off acoustic modes. 
This has implications as to its use in a variable area duct where modes may 
become cut off in the small area portion of the duct, or for situations where 
cut-off modes may be generated such as duct discontinuities or large velocity 
gradients. For single cut off mode propagation, Baumeister (ref. 81) found 
that the "steady state" impedance boundary condition produced acoustic reflec- 
tions during the initial transient which caused finite instabilities in the 
numerical calculations. In reference 81, extending the duct length to prevent 
transient reflections resolved this stability problem. 

Reference 81 also addresses the problem of hew long the transient calcula- 
tion must continue for the initial transient to die out in order to obtain a 
"steady state" solution. In agreement with analytical predictions (ref. 86), 
numerical calculation showed that the time to reach a steady state solution 
will be 


t* > L*/c* (30) 

where c* is the axial com.ponent of the group velocity. Since this 
g 

group velocity becomes very small near the cut off frequency, the transient 
solution will required long calculation times to resolve modes near cut off. 
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In proolems with duct orea variations, the finite element method is very 
convenient for "steady state" solutions of the wave equation. Unfcrtunately, 
the transient finite element method is inappropriate for two-dimensional proo- 
lems because the solutions are implicit requiring storage of the usual large 
matrix for each time iteration. Therefore, the transient finite difference 
theory must be adapted to duct problems with area variations. 

With area variations, the boundary can be located in an approximate manner 
with a uniform grid or can be located exactly by use of a variable spaced 
grid. Both approaches are inaccurate and are cumbersome to use. In reference 
82. White has developed a numerical mapping procedure which transforms a com- 
plex duct geometry into a simple rectangular form. Using the results of the 
numerical mapping, a transformed wave equation is solveo in the rectangular 
system with a standard urtifora rectaii 9 uidr 9 rid« 

Prof. White ot the University of Tennessee (Knoxville) under NASA Grant 
NAG3-18. is extending the procedure of reference 82 to the proolem of sound 
propagation from ducts into the far field. Figure 7 shows mapping for a 
straight two-dimensional duct into the far field while figure 8 illustrates 
the mapping for a variable area duct coupled to the far fiela. In both cases, 
the wave equation is solved in a rectangular grid system and transformed back 
to the appropriate grid locations shown in figures 7 and 8. 

In reference 83. as a check on both the steady state finite element and 
transient finite difference tneories, plane wave sound propagation was studied 
experimentdl ly in a rectangular duct with a converging-diverging area varia- 
tion ror no mean flow. The 0.5 area contraction was of sufficient magnitude 
to produce large reflections and induce some modal scattering. Figure 8 from 
reference 83 shows a comparison of White's (ref. 82) transient finite oiffer- 
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ence program and Astley-Eversman finite element program (ref. 74) with experi- 
mental data. Both theories accurately predict the standing wave pattern 
upstream of the test section (x < 0) and the pressure distribution through the 
area variation (0 £ x < 1). 

The transient theory has recently been applied to study the behavior of a 
sound pulse propagating through the shear layer of an axisymmetric jet 
(ref. 84). Both experimentation and numerical analysis show that in the low 
and intermediate frequencies the far field acoustic power exhibits a marked 
ampl ification as the flow velocity increases. The amplification is traced to 
shear noise terms which trigger the instability waves that are inherent within 
the flow. The experimental results were qualitatively in agreement with the 
numerical simulation. 

As discussed earlier in conjunction with the exit impedance, continuing 
the grid structure from inside the duct into the far field would simulate the 
actua*l dynamic process of reflection and transmission occurring at a duct 
lip. This is one of the strong motivations for Prof. White's study (ref. 82) 
anc the far field mapping shown in figures 7 and 8. In the far field, care 
must also be exercised to prevent false reflection generated at the termina- 
tion of the far field. Maestrello, Bayliss and Turkel (ref. 84) have includeo 
in their paper an excellent discussion of the problem involved in the far 
field termination as well as a correction to the exit i..pedance termination to 
improve the accuracy of the numeric:. I results. In addition, reference 85 con- 
tains a comprehensive literature summary of recent work on the external radia- 
tion boundary condition. 
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Status 


Although no documented computer program is available at the present time, 
the transient technique is relatively easy to program using tne publisheo 
theoretical results. The theory has been applied and checked for cases in- 
volving soft walls (ref. 61). uniform flow (ref. 62). variable area propaga- 
tion (ref. 79) and sheared mean flow (refs. 62 and 85). As yet, however, tne 
theory has not been checked in the combined case of a sheared mean flow with 
soft walls in a variable area duct. 

CONCLUDING REMARKS 

Tne finite difference and finite element theories are iaeally suiteo tor 
predicting sound propagation in ducts particularly in low frequency applica- 
tions such as mufflers, expansion chambers, and exhaust ducts of turbofan 
engines. Using available computer programs, souna attenuations can now be 
easily and precisely calculated for ducts with a variety of complexities, suen 
as variation in the wall liner impedances, axial area changes or large varia- 
tions in the mean flow field. Of course, the sound source distribution and 
acoustic liner parameters must be accurately known. 

On inn thnorntlcal sloe, et tne present tin*, reseeren priorities tnclooe 
extenoin, tne numerioel tneories to nigner freguencies, coupling tne Internal 
portion or tne ouct to tne far flelo. ano incluoing non-linear effects. 
Experlwntolly. testing eitn large area ano luean floe variations Is continuing. 
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Figure U • lumped-parameter representation d continuous acoustic flow field. 
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